Propagation Of Waves In Periodic-Heterogeneous Bistable Systems 
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Wave propagation in one-dimensional heterogeneous bistable media is studied using the Schlogl model as 
a representative example. Starting from the analytically known traveling wave solution for the homogeneous 
medium, infinitely extended, spatially periodic variations in kinetic parameters as the excitation threshold, for 
example, are taken into account perturbatively. Two different multiple scale perturbation methods are applied 
to derive a differential equation for the position of the front under perturbations. This equation allows the 
computation of a time independent average velocity, depending on the spatial period length and the amplitude 
of the heterogeneities. The projection method reveals to be applicable in the range of inteiTnediate and large 
period lengths but fails when the spatial period becomes smaller than the front width. Then, a second order 
averaging method must be applied. These analytical results are capable to predict propagation failure, velocity 
overshoot, and the asymptotic value for the front velocity in the limit of large period lengths in qualitative, often 
quantitative agreement with the results of numerical simulations of the underlying reaction-diffusion equation. 
Very good agreement between numerical and analytical results has been obtained for waves propagating through 
a medium with periodically varied excitation threshold. 
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I. INTRODUCTION 



A. Heterogeneities In Reaction-Diffusion Systems 



Spatio-temporal patterns of Reaction-Diffusion Systems 
(RDS) are of fundamental interest in many chemical [1] and 
biological systems |2]. Prominent examples are the Belousov- 
Zhabotinsky reaction (BZR), action potential propagation in 
cardiac tissue and chemical catalysis. A large variety of pat- 
terns can be found, e.g. traveling fronts and pulses in one 
spatial dimension, spiral waves and traveling spots in two and 
three spatial dimensions yl 01 ■ Most of the models describing 
these effects are assumed spatially homogeneous, although at 
least biological systems are intrinsically heterogeneous. Re- 
cently there is an increasing interest in heterogeneous RDS 
where the diffusion coefficients, the reaction rate constants or 
other important parameters as the excitation threshold depend 
explicitly on space. Effects of heterogeneities on traveling 
wave solutions of RDS reach from reflection and diffraction 
to velocity overshoots (the velocity of the wave is larger in the 
case of heterogeneities than without them) and propagation 
failure or oscillatory pinning [5-101. In two spatial dimen- 
sions, breakup of plane waves into spiral waves and spatio- 
temporal chaos can occur Illll[l2ll . 

The BZ reaction is an example for an experimental realiza- 
tion of a heterogeneous RDS. The use of masks and litho- 
graphic techniques permits the introduction in the reaction of 
patterned illumination and patterned distribution of catalyst, 
respectively. A spatially varied intensity of applied light cor- 
responds to a spatial variation of the excitation threshold of the 
system. A mathematical model for this reaction is the modi- 
fied Oregonator model 1 13]. The velocity of pulse propagation 
in one spatial dimension under the influence of a spatially pe- 
riodically rectangularly varied excitation threshold was stud- 
ied numerically in f 14] and revealed a velocity overshoot at 
small period lengths. An analytical investigation by Keener 
^ based on the averaging theorem for the Schlogl model 
with a spatial variation of a reaction rate in form of a Dirac 



comb showed a large velocity overshoot for period lengths of 
the heterogeneities smaller than the frontwidth. We derive a 
slightly generalized version of Keener's method in section III. 
An alternative perturbation method was applied to scalar and 
multicomponent RDS by Bode JT], [^ and also by Nishiura 
la uM with mostly bump-type heterogeneities. In section II 
a different form of this perturbation method, called projec- 
tion method throughout this publication, is proposed, which 
allows the investigation of the effect of heterogeneities with 
medium and large period lengths on front propagation. In the 
form we use this method it was also applied in e.g. IUSI - UTII 
to traveling fronts in stochastic bistable media. Both meth- 
ods are applied to various variations of kinetic parameters of 
the Schlogl model in section IV, and compared with numerical 
results in section V. 



B. The Schlogl model 

This scalar RDS, proposed by Zeldovich-Frank- 
Kamenetsky [18] as a model for front propagation, and 
later applied by Schlogl [19] as a model for a non-equilibiium 
phase transition, also known as bistable equation, is given in 
its general form as 



d,u = Dd^u + R{u) 



(1) 



with a reaction function 



R{u) = —k{u — Ml) (m — M2) (m — M3) , < Ml < M2 < M3- (2) 

The parameters mi, M3 are stable fixed points and the param- 
eter M2 is an unstable fixed point which corresponds to the 
excitation threshold of the system and D is the diffusion co- 
efficient, k is called reaction coefficient and has a dimension 
[time]^ . It is a measure of the intrinsic time scale on which 
the reaction takes place. The traveling front solution of (|2]) is 



- Mi+M3 + (Mi-M3)tanh -W — (m3-mi)^ J J (3) 



with ^ —x — ct and a velocity 



Dk . ^ . 



(4) 



The front width / of the traveling wave solution can be defined 
aslfiil 



(5) 



4V2D 



\/k{uj — Ml) 



For every choice of the value of the excitation threshold U2 
the front has a certain velocity c ~ 0(112) but the same front 
profile Uc; which shows no explicit dependence on U2- This is 
a peculiarity of the Schlogl model. The general Schlogl model 
(|2]l can be cast, without loss of generality, into the form 



dfU — d^ u ~ u (u — U2) (m — 1) , < M2 < 1. 
The traveling wave solution, with boundary conditions 

limM(jc,f) = 0, lim u(x,t) — l, 

simplifies to 

/ 1 



Uci^)-- 
with a velocity 



1 - tanh 



V2V2' 



1 



A. 

,s/2 



c^-^{l-2u2). 



(6) 



(7) 



(8) 



(9) 



All analytical computations are done with the simpler form 
(|6]l of the Schlogl model. 



C. Harmonic mean velocity 

The appropriate average speed for a wave traveling with a 
space dependent velocity is the harmonic mean of the velocity. 
Suppose a wave travels a distance L/2 with velocity ci and the 
same distance L/2 with velocity 02- The total time it takes the 
wave to travel through both regions is 



L / 1 1 

T = Ti+T2 = -{- + - 
2 \ci C2 

and the harmonic mean velocity is 



L 2 

T l/ci + l/c2 



(10) 



(11) 



If the velocity is assumed to depend on the space coordinate x 
in a periodic way with period length \,c — c{x), one can ap- 
proximate the average velocity over an arbitrary period length 
L by dividing Linn pieces 



^har> 



1:1^,1 1 c{xi/L) 



(12) 



Consider a traveUng wave with speed c depending on a pa- 
rameter s,c = c{s). If this parameter is spatially varied in the 
form of an infinitely extended periodic function with period L, 
s = s [x/L) , the harmonic mean velocity can be approximated 

as 

The underlying assumption is that the front instantaneously 
adapts its velocity when transiting from one region of space 
to the other, distinguished by the different values s [xi/L] of 
the parameter s. In the limit of infinitely many pieces the sum 
becomes an integral and the harmonic mean velocity is 



Che 



= 1/ 



1 



c{s{x)) 



Ax 



(14) 



For the more realistic case that it takes the front some tran- 
sient time to adapt its velocity when transiting from one region 
to the other, one can state that c/,a„ij (O still gives a good 
approximation for the average speed, if the transient time is 
small compared to the traveling time through one period of 
the spatial heterogeneities. This is always the case for the 
limit of an infinite period length of the spatial heterogeneities. 
Thus one expects the harmonic mean velocity (fT4l l to give the 
approximate average velocity for a wave traveling through a 
periodic medium with large period lengths. 



II. PROJECTION METHOD 

Consider an unperturbed scalar RDS in one spatial dimen- 
sion, Eq. ([T]i (with D=\) and a traveling wave solution Uc {^ ) 
with constant velocity c. A perturbation K {u,x,t) , depending 
on space and time as well as on m, is introduced. K is multi- 
plied by e, which serves as the small parameter for the pertur- 
bation expansion and is set to e = 1 at the end of the compu- 
tations. The perturbed RDS, written in the co-moving frame 
E, =x — ct of the unperturbed RDS dU and with m = m (^ , f ), is 

dtu = dpU + cd^u+R{u) + eK{u,S, +ct,t). (15) 

The introduction of an additional heterogeneous diffusion co- 
efficient is straightforward [21], but not done here. Under the 
perturbation K, the approximate solution of the perturbed RDS 
(fTSl l is assumed to be 



u{^,t)^Uci^) + Eui^,t) 



(16) 



Inserting (fTSI l into the unperturbed RDS ([T) and expanding R 
in powers of e leads, in order e' , to 

d,u^d^a + cd^u+R'{Uci^))a^^u. (17) 



The operator ^ is a linear differential operator with eigen- 
values A which determine the stability of the traveling wave 
solution Uc- Under the condition of a translationally invariant 
reaction function in ([T]i one can prove the existence of a cer- 
tain eigenfunction S (^ , f ) = [// ((^ ) of ^ corresponding to the 
eigenvalue Aq = which is called the Goldstone mode. If the 
wave is stable, this is the largest eigenvalue. The function u 
can be split up in a part parallel to the Goldstone mode and 
and a part orthogonal to it (with /?, q arbitrary constants), 



u{^,t)=pU'Ai)+qv{^,t) 



(18) 



with 



{W\^),v{^,t)) = jw\^)v{^,t)A^=Q. (19) 



Eq. ( fT9] l is called projection condition, where 



^„f)=lg{^)f{^)d^ (20) 



is the inner product in function space and W^ is the Goldstone 
mode of the adjoint operator of J^, J^'^W'^ = 0, and will be re- 
ferred to as the adjoint Goldstone mode. The Goldstone mode 
U'^ leads to a small shift of the traveling wave solution, 

f/,(<^) + pf/;(^) + <?v(^,r)«C/,(<^+p) + <?v(^,f), (21) 

so the projection condition ( fT9] l can be interpreted as saying 
that V does not take part in a shift of the traveling wave solu- 
tion, it only leads to a deformation of the wave profile. 
A slow time scale T = et is introduced and the time derivative 
transformed accordingly 



d,^d, + dfTdr = d, + edj. 



(22) 



From now on, T and t are treated as independent of each other, 
as it is the usual procedure in multiple scale perturbation the- 
ory. The ansatz for the solution of the perturbed problem (ITsT l 

is 

M(^,f,r) = f/,(^+p(r)) + ev(^,r,r), (23) 

together with the projection condition 

(w'(<^+/,(r)),v(<^,r,r)) = 

jwH^+p{T))v{^,t,T)A^^O, (24) 

where the correction to the position of the front under per- 
turbation p (r) is constant on the original time scale t but 
depending on the 2nd timescale T. The functions R{u) and 
k{u,^ +ct,t) are expanded in powers of £. In order e' one 
gets 

-d,v{^,t,T)+^vi^,t,T)^ 
-K{Uc{^+p{T)),^+ct,t)+p'{T)U^{^+p{T)), (25) 



a linear PDE for the correction v of the front shape under per- 
turbation with an inhomogeneity on the rh.s. Eq. (l25T l is pro- 
jected onto the adjoint Goldstone mode, and by applying a 
variant of the projection condition (|24] | one can eliminate one 
term to get 

{wH^+piT)),^v{^,t,T)) = 

-{wH^+p{T)),KiUA^+p{T)),^+ct,t) 

-p'{T)U',{^+p{T))) (26) 

This turns out to be a solvability condition (or Fredholm alter- 
native, see [22]) for v, which guarantees a bounded solution 
for V only if the rh.s. of ( |26] | is zero. With the adjoint Gold- 
stone mode of a scalar traveling wave solution Uc (<^) in one 
spatial dimension. 



W'(<^)=e^-^f/.'(^) 



(27) 



and a coordinate change ^ -^ ^ + p (T) , the solvability con- 
dition reads as 



je'^U',{^){K{U,.{^),^+ct-p{T),t) 



-p'{T)U^{^))d^=Q. (28) 

This is a differential equation for p{T) , the correction to the 
position of the front. Rescaling p (T) to the original time t by 
introducing a new function 



^it)^ct-p{T) 



and using ( |22] ) gives 
d 



^H')-c-e-piT). 



(29) 



(30) 



Thus one has derived the ODE for the position of the traveling 
wave under perturbation 



-0(r)-c + e0f((/.(f)) 



with 



Kc^ je^^Ul^i^)^^. 



(31) 



(32) 



The advantage of introducing a new function , (|29] l, is that if 
the perturbation K has no explicit time dependence, one gets 
an autonomous ODE for (f ) which is usually much easier to 
solve. If the rh.s. of dSTl i is zero, propagation fails. 
In the case of a time independent perturbation K with period 
length L 



k{u,x) = k(u,x + L) , 



(33) 



one can compute a time-independent average velocity Cavg 
over one period of the heterogeneities. The travehng time T^ 
is the time it takes the front to travel through one period of the 
spatially periodic heterogeneities 



Tc 



d(j) 



= L 



C + £0[(0) J C + £0f(0L) 





d0 



(34) 



and the average velocity is 

Cavg = — = 1/ 



d^ 



c + e©f(0L)' 



(35) 



III. AVERAGING METHOD 

Keener developed a method based on the averaging theorem 
II22I I23I1 to derive an ODE for the position of the front of a 
RDS with a heterogeneous reaction function of the form 



R{u,x) = I 1 



f (u) — a{u) , 



(36) 



where / is an arbitrary nonlinear function |J5ll- He also treated 
the case of a heterogeneous diffusion coefficient 12411 with a 
similar method. The spatial heterogeneities have to be peri- 
odic and the period length L is used as the small parameter for 
the perturbation expansion. 

Written as a system of two differential equations with reaction 
function ( [36] ), ([B is 



dxU — v, 

5,y==a,M-(l+/(0)/(M)+fl(M). 



(37) 



In contrast to Keener's approach pi], which assumes a linear 
dependence of a on m, here the function a (w) is allowed to be 
an arbitrary, possibly nonlinear function. The function g' (x) 
denotes the heterogeneities with period one and zero mean: 

1 

{g'{x))^Jg'ix)dx^O. (38) 



A second length scale T = x/L is introduced and the space 
derivative transformed accordingly, dx -^ dx + j"?!- An exact 
change of variables from m, v to new variables y, z is applied 
to (|J7b to get, after expanding nonlinear terms in (l37t and 
the unknown functions y, z in powers of L, a homogeneous 
averaged system 

dtyQ{x,t)^djyQ{x,t)+f{yQ{x,t))-a{yQ{x,t)) (39) 

in order lP and linear inhomogeneous PDEs for the new vari- 
ables in higher orders of L. Eq. (l39T l must have an analytically 
known traveling wave solution yo [x^t) — Uc {x — ct) . Trans- 
forming into the co-moving frame ^ = x — (p (t) with an un- 
known, time dependent position of the front ^ (t) , and apply- 
ing a solvability condition to the linear PDE obtained in the 
next non- vanishing order of L yields an ODE for (j) {t) . For 
more details, see Keener's publication l^. 



A. Averaging in 1st order 

The exact change of variables is 

u{x,t,T) =y{x,t,x) 
v{x,t,T)=z{x,t,T)-Lg{T)f{u), (40) 

where g is the anti derivative of g\ The ODE for the position 
of the front under perturbation is 



dt 



:c + L0f(0(r)) 



(41) 



with 



and Kc given as above (|32] |. 



B. Averaging in 2nd order 

For second order averaging a different exact change of co- 
ordinates is applied which does not only eliminate all het- 
erogeneous terms in order LP , but also all terms in order L' 
and yields a linear inhomogeneous PDE for the new variables 
y2 , Z2 in order L^ . 
We follow Keener and use as exact change of coordinates 

u{x,t,T)=y{x,t)~L^G{T)f{y), 
v{x,t,T)=z{x,t)-Lg{T)f{y)+L^G{x)f'{y)yx 



+ L'giT)GiT)f'iy)fiy), 
where G is the anti derivative of g. If 

lim /(C/,(^))=0, 



(43) 



(44) 



the integration constants for g and G can be chosen so that the 
mean values {g{x)) and {G{x)) are zero. If not, one has to 
choose 



{G{x))^{g^{x))lim - 



f'{Uc) 



i^-^f'iUc)-a'{Uc) 



(45) 



The ODE for the position of the front under perturbation in 
2nd order averaging is 



#(0 

dt 



:C + L"©f(0) 



(46) 



with 



0^(0) = -yJs' l^^^)riuc)fmuy^d^ 



^'^)^(^^'^(^(^^)^''))'^^-^^^^ 



The explicit dependence on the homogeneous part a [u) of the 
reaction function can be eliminated in both cases W\\ . ( |46] l. 
which is also the reason why these results are the same as de- 
rived by Keener [5] for a linear function a (u) . Eq. ( |42] | and 
(|47] | depend implicitly on a through the traveling wave solu- 
tion of the homogeneous case Uc {£, ) ■ 

To compute an average velocity Cavg over one period length of 
the heterogeneities, one proceeds analogously to the projec- 
tion method (l35b. and derives 



^avg 







d(j) 



-L^&ji^L)' 



(48) 



C. Equivalence of 1st order averaging and projection metliod 

With one partial integration and assuming that 

g ( ^-j^ J / (Uc) Ul-e^^ vanishes at the boundaries, the ODE 

for the position of the front derived in first order averaging 
dTTT i becomes 



d0(O 
dt 



K, 



f{Uc)Uy^d^ (49) 



If the periodic-heterogeneous reaction function can be split up 
in the form we assumed for the averaging method, R {u,x) ~ 
f {u) - a{u) + g' {x/L) f {u) =/?(m) + k-(m,.^:), then the ODE 
for the position of the front derived with the projection method 
(ISTT i becomes 






yJ g'[^^)fmuy^A^- (50) 



For evaluation of (ISOl l, e is set to e = 1 and one realizes that 
dSOl l and ( |49] l are equal. Note the differences in the approaches 
of these two multiple scale perturbation expansions: for the 
projection method, a second time scale T = ef is introduced. 
For the averaging method a second space scale T = x/L is 
introduced and the heterogeneities are restricted to periodic 
ones with small period lengths. Thus it can be shown that 
the validity of equation dTTT l for the position of the front de- 
rived in 1st order averaging can be extended to arbitrary large 
period lengths and even to non-periodic heterogeneities. In 
fact, in the case of the Schlogl model, it fails for small period 
lengths but gives good results for large period lengths, as will 
be shown later. 



IV. ANALYTICAL RESULTS FOR THE SCHLOGL 
MODEL 

A. Projection metliod for a variation of «i , m2 , "3 and k 

An infinitely extended sinusoidal variation of the excitation 
threshold U2 of the Schlogl model is considered. 



M2 {x) = U2 + A sin {2nx/L) 



(51) 



where A is the amplitude of the spatial variation and L its pe- 
riod length. The heterogeneous reaction function is 

R{u,x) = —u{u— {u2+Asm{2nx/L))) (m— 1) 

= — m(m — M2) {u— 1)+Asm{27rx/L)u{u— 1) 
^R{u) + k{u,x). (52) 

Similarly, one introduces the heterogeneities for the variations 
of Ml , M3 and k. The derivation of the ODE for the position of 
the front under perturbation can be done simultaneously for 
all four variations by introducing a general perturbation 

k{u,x) = (m — Zi) (m — Z2) {Z4U+Zj)Asm{27tx/L) , (53) 

where one has for the different variations 

1. Ml (x) =Asm{27tx/L) : 

Zl = 1 , Z2 = M2 , Z3 = 1 , Z4 — 0, 

2. M2 (x) — M2 + A sin {2nx/L) : 

Zi = l,Z2-0,Z3 = l,Z4 = 0, 

3. M3 (x) = 1 +Asin(27rx/L) : 

Zi = 0, Z2 = M2, Z3 = 1, Z4 = 0, 

4. k(x) = l+Asm{2nx/L) : 

•Zl = 1 , Z2 = M2 , Z3 = 0, Z4 = — 1 . 

Note that the amplitude A is restricted to values for which the 
condition < mi < M2 < M3 and A: > is fulfilled locally, e.g. 
for a variation of M2 : mi < M2 (x) < M3 for all x. This implies 
that the computation for a variation of mi cannot be done with 
the form (|6]l of the Schlogl model, where mi = 0, but the gen- 
eral form (|2]i has to be used instead. 

The differential equation for the position of the front (f) , see 
(EB, is 



X (f/, (^) -Zi) (f/,(^) -Z2) (Z4f/, (<^) +Z3)d^ 

= c + sCi(c2sin(^)+C3Cos(^ 

(54) 

The values of the constants Ci , C2 and C3 are given in the ap- 
pendix, see ( l76b . The average velocity Cavg can be computed 
according to the formula (l35T l to get 



(55) 



C.v, = V^2_£2c2(c2 + C32). 



The ODE for the position of the front (l54l with the initial 
condition (f)(0) —0 can be solved to give 



C. 2nd order averaging for a sinusoidally varied excitation 
tlireshold 



^,, L / c + eCiC 

mlt) — — ai-ctan 7- r ^ 

TT \cOt{lCa,gt)ca,g-eCiC2 



(56) 



which describes the position of the front for one spatial period 

-§<0(r)<§. 



of the heterogeneities, — 4 < (f) < ^ 



Propagation failure occurs if -j-^ (f ) = 0. With the help of the 
r.h.s. of ( |54] | one derives the relation ll25ll 



^avg ■ 







(57) 



as the condition for propagation failure. 
Computation of the ODE for the position of the front d3TT i can 
be done for periodic arbitrary shaped heterogeneities by ex- 
panding k{u,x) in a Fourier series in x. This was done for 
variations of all kinetic parameters in form of an infinitely ex- 
tended rectangular function, but results are not shown. 



B. Failure of projection method for small period lengths 

As was already mentioned by Keener 0. 12411 . 1st order av- 
eraging fails for small period lengths and smooth nonlinear 
functions f{u) . This can be shown for the Schlogl model by 
estimating the dependence of the coefficients C1C2 and C1C3 
in (|54| on L, see ( |76t . For small L, one can approximate 



cosh 



sinh 



a 






(58) 
(59) 



in C2 , C3 , and in the denominator of Ci 

2 



1 (iVlcn) 



cosh 



4V2n- 



< 



1 + cosh 



'4^271:2 



(60) 



1 4V1k^ 



By keeping only the largest terms up to order L' in Hi, H2 
(dTTl), dTSll), one derives 



The heterogeneous reaction function is cast in the form nec- 
essary for averaging (l36t 

R{u,x) = —u{u— {u2+Asm{27tx/L))) (m— 1) 



IH sm{27tx/L) ] u(u~ I) U2 — U (m— 1) 

M2 



= 1 



+ 8'{l))f{u)-aiu). 



(62) 



To derive the differential equation (|46t one has to solve the 
integrals arising in ©2 , ( |47T i. The integration constants for the 
anti derivative g of g' and G of g can be chosen so that g and 
G have vanishing mean, (|44] |. The function 02 (0) can be 
computed analytically to get 



@2i<P) =Hism 



m 



-H2Cos(^) 



-//3sm(^— j 



//a cos 



m 



H5. (63) 



The values of the constants Hi up to H^ can be found in the 
appendix, see dSOl l. The average velocity (|48T l is computed 
numerically, because analytical computation is possible, but 
very tedious. 

Computation of the ODE for the position of the front (ISTT i can 
be done for periodic arbitrary shaped variations of the excita- 
tion threshold U2 by expanding g and G in ( |47] i in a Fourier 
series in x. This was done for a rectangular variation of the 
excitation threshold, but results are not shown. 



D. 2nd order averaging for a sinusoidally varied fixed point 
parameter 1/3 



The heterogeneous reaction function (l36t is 

R{u,x) = —u{u — U2) (m — (1 + A sin {2nx/L))) 

= (1 +Asm{27tx/L))u{u — U2) — m (m — M2) 
'l+fi'('^\)f{u)-a{u). (64) 



2 2 
c — c 



e^C\{Cl + Cl) 



A\/2lt^ 



256A^e'^e t~ %^ sin 



U2c%\ 



47:^sm^(V2c7:]zl 



{c-2c^fL^ 

2;rcos ( \/2ck] Z4 + Lsin I \f2cn 

X ( (4c + V2 (2Zi + 2Z2 - 3)) Z4 - 2 V2Z3) ) J + O (L^) 

(61) 

The term e-4\/27r /l vanishes faster than any polynomial order 
of L, and is called transcendentally small by Keener 



The integration constant for G is determined so that (G) = 

8^43»^-5) ' 'S3'- "^^^ ^^^ ^°'' *^ position of the front has the 
same form as for a variation of M2, 



0^W=//,sin(^)+i/2Cos(^) 

+ H^%m f ^ j +//4COS y-^ ) +//5, (65) 



The constants H\ up to H^ are listed in the appendix, see 
The average velocity ( |48] | is computed numerically. 
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Figure 1 : Ratio of velocities Cavg/c plotted over tlie ratio of period 
length L and frontwidtli I for a sinusoidal variation of excitation 
threshold M2 W = 0.35 +Asin(2;rx/L) . Projection method (green 
solid line) shows excellent agreement with numerical results (red 
crosses), 2nd order averaging (blue dashed line) fails for large pe- 
riod lengths. The harmonic mean velocity (purple dash-dotted line) 
gives the approximate average velocity for large period lengths. 



Figure 2: Velocity overshoot for a sinusoidal variation of excitation 
threshold U2 (x) = 0.35 +Asin(2;rx/L) , predicted by 2nd order av- 
eraging (blue dashed line) in qualitative agreement with numerical 
results (red crosses). Projection method (green solid line) fails for 
small period lengths. 



V. COMPARISON OF ANALYTICAL AND NUMERICAL 
RESULTS 

The PDE for the heterogeneous Schlogl model is solved 
numerically with a simple finite differences Euler forward 
scheme and the average velocity is obtained by a linear fit 
of the space over time data over an integer multiple of period 
lengths. Data near to the boundaries have to be neglected be- 
cause of boundary effects. 



small period lengths but fails for large period lengths, as could 
be expected because the period length is used as the small pa- 
rameter for the perturbation expansion (although 1st order av- 
eraging does not fail for large period lengths). 



A. Variation of excitation threshold M2 

A sinusoidally varied excitation threshold ui {x) — U2 
Asm (^) leads to an average velocity 



4A- 



((m2- 1)2^2 



L-' + lK' 



(2U2 — 3U2 + U2) L^ 



{uIl^ + iTZ^) (l2 (1 - 2u2f + inA sin^ (2u2n) \ ' 



cosh fi^")- cos (4M2?r) 



(66) 



obtained by the projection method according to the formula 
dSSl l. For intermediate and large period lengths, the agree- 
ment between ( |66] l and numerical results is excellent, see Fig. 
[T] The averaging method in 2nd order gives good results for 



B. Velocity overshoot 



For period lengths smaller than the front width, the nu- 
merical results for the sinusoidally varied excitation threshold 
show a small velocity overshoot. The results obtained in 2nd 
order averaging predict this overshoot qualitatively at the right 
size of period lengths. Eq. (l66T l shows a plateau-like behavior 
indicating the transcendental smallness of the results obtained 
with the projection method for small period lengths, see Fig. 



For a sinusoidal variation of the fixed point parameter 
M3 (x) = 1 + A sin ( 2p ) , a large velocity overshoot is found 
numerically, again predicted qualitatively by 2nd order aver- 
aging and missed by the projection method, which gives an 
analytical solution for the average velocity 
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Figure 3: Velocity overshoot for a sinusoidal variation of 1(3 (x) = 
1 +Asm{27tx/L) and 112 = 0.35 predicted by 2nd order averaging 
(blue dashed line) in agreement with numerical results (red crosses). 
Propagation failure is predicted qualitatively by projection method 
(green solid line), which fails for small period lengths. 



Figure 4: Propagation failure for a sinusoidal variation of excitation 
threshold M2 (x) = 0.35 +A sin {2nx/L) occurs if the average velocity 
obtained by numerical simulation (red crosses) drops to zero and can 
be predicted by projection method (green solid line) and qualitatively 
by 2nd order averaging (blue dashed line). 



(L^{l-2u2f + SnA {uJL^ +27:^) sin^ {2u27:)\ ' 



cosh ( ^4^)- cos (4M2?r) 



(67) 



large enough, propagation failure occurs for all period lengths 
larger than a certain critical period length, see Fig. |4] For the 
choice of parameters shown in Fig. |4] the result of 2nd order 
averaging can predict the propagation failure, but if propaga- 
tion failure occurs at larger period lengths, it fails badly to 
do so. Comparison of ( |68] | with numerical results shows very 
good agreement, see Fig. |5] 

For a sinusoidal variation of the reaction coefficient k (x) — 
1 + Asin (^) a different behavior occurs. The average ve- 
locity obtained with the projection method is given as 



C. Propagation failure 



With the condition for propagation failure, Q,vg = 0, the 
results for the average velocity obtained with the projection 
method allows to determine the critical amplitude for which 
propagation failure occurs. For the case of a sinusoidally var- 
ied excitation threshold, this gives 



l3 , /cosh f ^4^ ) - cos (4;rM2) 



Acrit.U2 (.^j — 



2V2VL2(m2-1) +27r2 



CSc(27rM2)(l-2M2) (1-M2)m2 ,,„. 



L^uI + 2k^Jl^ (1 - 2m2)^ + 8;r2 



which is a monotonically decreasing function for all values of 
the system parameter U2 ■ This means that for a fixed amplitude 
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4m^(m2(2m2-3) + 1) 
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cosh(^4^)-cos(4M27r; 



(69) 



and shows, for a certain range of the system parameter U2 and 
in qualitative agreement with numerical simulations, a min- 
imum for a finite period length, see Fig. |6] For a value of 
the amplitude large enough and slightly different values of the 
system parameter U2, the projection method predicts propa- 
gation failure occurring for an interval of peiiod lengths, see 
Fig. |7] and propagation is possible for larger and smaller pe- 
riod lengths. Comparison with numerical results shows quali- 
tative agreement, see Fig. |8] The critical amplitude for which 
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Figure 5: The critical amplitude for which propagation failure oc- 
curs for a sinusoidal variation of the excitation threshold U2 (x) = 
U2 -\-A sin {Inx/L) , obtained by the projection method, is a monoton- 
ically decreasing function for all values of 112 and shown for U2 = 0.3 
(blue dashed line) and «2 = 0.4 (red dash-dotted line). The results 
for a value of U2 = 0.35 (green solid line) are compared to numerical 
simulations (red dots). 



Figure 6: Ratio of velocities Cavg/c plotted over the ratio of period 
length L and front width / for a sinusoidal variation of k{x) = 1 + 
A sin [iTtx/L) and U2 = 0.38. The average velocity shows a minimum 
and the numerical solution (red crosses) approaches the harmonic 
mean velocity (black dashed line) from below. The solution obtained 
by projection method (green solid line) approaches a different limit 
(purple dot-dashed line). 



propagation failure occurs is given as 
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Acrit.k (,-^j 



cosh(^ 



COS (4-71112) 



L^{u2-lf + 2TZ^ 



csc(27rM2)(l -2m2) (1-M2)m2 



(70) 



and is shown for different values of U2 in Fig. |9l There one 
can see again the fact that propagation fails for an interval 
of period lengths: the critical amplitude has a minimum for 
a finite period length and is not a monotonically decreasing 
function. 

For most cases of a variation of M3, the critical amplitude 
for which propagation failure occurs is a monotonically de- 
creasing function, but for a very small range of the system 
parameter M2j propagation failure can occur for an interval of 
period lengths (not shown). 



D. Front velocity in the limit of large period lengths 

The average velocity obtained with the projection method 
allows to determine the limit for large period lengths. For the 
case of the sinusoidal variation of excitation threshold M2, this 
limit is given as 



lim c, 



uvg, U2 



= Vc2-2A2, 



(71) 



which agrees with the harmonic mean of the velocities com- 
puted according to (fT4l i. The numerical solution approaches 
the limit (ItT] ) from above, as shown in Fig. [T] The limit of the 
average velocity 



lim Cavg.u-i 



- J A 



2^ ' 



(72) 



for the sinusoidal variation of M3 also agrees with the harmonic 
mean of the velocities and with numerical results (not shown). 
For the sinusoidal variation of reaction coefficient k the 
limit for large period lengths of (|69^ is 



lim c, 



avg.k ■ 



-^T. 



(73) 



which does not agree with the harmonic mean of the veloci- 
ties. 
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(74) 



where K (x) is the complete elliptic integral of the first kind. 
The numerical solution approaches the harmonic mean of the 
velocities ( f74b from below, see Fig. |6l Fig. |8] The velocity 
of the unperturbed general Schlogl model has the same square 
root dependence on the reaction coefficient k as on the diffu- 
sion coefficient D, see dUl. It follows that the limit of large pe- 
riod lengths of the average velocity for the case of a sinusoidal 
variation of the diffusion coefficient D{x) = \+A sin {iKx/L) 
is the same as (|74] i, which was checked by numerical simula- 
tions (not shown). 
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Figure 7: The ratio of velocities Carg/c either shows an interval of 
period lengths for which propagation failure occurs or a minimum 
for a sinusoidal variation of k{x) = i+A sin {iTtx/L) , depending on 
the value of U2 = 0.35 (blue dashed line), U2 = 0.425 (red dot-dashed 
line), U2 = 0.38 (green solid line). 



VI. DISCUSSION 

Infinitely extended spatially varied reaction parameters of 
the Schlogl model are considered and the effects on the prop- 
agation velocity are studied. The applied perturbation meth- 
ods seem to work best for a variation of excitation thresh- 
old U2, the reason is probably that the front profile ([3) does 
not depend on this parameter. The projection method works 
worse for a variation of the reaction coefficient k than for all 
other variations, even not predicting the correct limit of the 
average velocity for large period lengths. This could be con- 
nected to the fact that the velocity of the homogeneous case 
(|4]i shows a linear dependence on the fixed point parameters 
but a square root dependence on k. The ODE for the position 
of the front obtained in 1st order averaging is equivalent to the 
one obtained with the projection method, and both fail gener- 
ally for small period lengths due to the transcendentally small 
dependence of ©i (0) on the period length. This causes the 
plateau in the plots of the average velocity for small period 
lengths in all solutions obtained with the projection method. 
The solutions obtained in 2nd order averaging agree qualita- 
tively with the numerical simulations and can predict the size 
and the value of the period lengths for which velocity over- 
shoots occur, but generally fail for large period lengths. A 
small velocity overshoot up to 1.5% is found for a variation of 
the excitation threshold U2 with period lengths slightly smaller 
than the front width. A similar velocity overshoot was found 
for period lengths for a variation of the excitation threshold 
in the modified Oregonator model 1 14). For a variation of M3, 
a larger velocity overshoot up to 25% is found, which occurs 
at period lengths of approximately the same size as the veloc- 



Figure 8: Interval of period lengths for which propagation fail- 
ure occurs for a sinusoidal variation of reaction coefficient k (x) = 
1 +As'm{2nx/L) , predicted by projection method (green solid line) 
in qualitative agreement with numerical results (red crosses). The 
limit for large period lengths of the numerical results is given by the 
harmonic mean of the velocities (black dashed line), the analytical 
solution approaches a different limit (purple dot-dashed line). 



ity overshoot in the case of the variation of U2- Propagation 
failure occurs in both cases for amplitudes large enough and 
all period lengths larger than a certain critical period length. 
For a variation of the reaction coefficient k, and for a small 
range of values of the system parameter U2 in the case of a 
variation of M3, we find an interval of period lengths for which 
propagation failure occurs. All computations were done for a 
sinusoidal as well as for a rectangular variation of the parame- 
ters, which show qualitatively the same effects. The amplitude 
A and the period length L are more important in affecting the 
front velocity than the shape of the heterogeneities. 



Appendix A 

The ODE for the position of the front for a sinusoidal vari- 
ation of Ml, M2, M3 and k obtained with the projection method 

is 



X {Uc (^ ) - Zi ) {Uc (^ ) - Z2) (Z4M + Z3) d^ 
:c + eCi(C2Sin( — p^j+C3Cos' ^^ ' 



with 



(75) 
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icrit.i: 




H2 ^ LSn^ (2V2Z3 -Uc + Vl {Til +2Z2-3)]Z4^ 
^L^7t((6V2{2Zi+2Z2-3)c^ + {24Zi{Z2-l)~24Z2)c 
+22c + Sc^ + V2 (Zi (4 - 6Z2) + 4Z2 -3))Z4 
-2 (6V2c^ + 12 (Zi +Z2 - l)c 
+\/2(-3Z2+Zi(6Z2-3)+2))z3). (78) 



Appendix B 

With the help of the averaging method in 2nd order an ODE 
for the position of the front for a sinusoidal variation of the 
excitation threshold is derived 



Figure 9: The critical amplitude for which propagation failure occurs 
obtained by projection method for a sinusoidal variation of k {x) = 
1 +A sin {27tx/L) shows a minimum for all values of 112 = 0.35 (blue 
dashed line), 112 = 0.425 (red dot-dashed line) and 112 = 0.38 (green 
solid line). 
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Appendix C 



For the sinusoidal variation of the fixed point parameter M3, 
2nd order averaging gives an ODE for the position of the front 



©f ((^) = Hi sin {^j +//2COS {^j 

+ H3sm(^]+H4Cos(^]+H5, (83) 
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and Di, D2 are given as in 
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determine the value of ^ for which the r.h.s. of ( 154) attains its 
minimum, substitute this value back into the r.h.s. of ( I54t and 
determine a relation between the parameters under which the 
r.h.s. of ( I54t is zero. 



